home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / TRY / my_fft.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  11.8 KB  |  534 lines
  1. #define MAX_SIZE 1111
  2. #define MAX_THREADS 8
  3. *****************************************************************************
  4. *                                        *
  5. *    1D (slow) Fourier Transform                        *
  6. *                                        *
  7. *****************************************************************************
  8.       subroutine dft1du( plusminus, n, a, inc)
  9.       double precision a(*)
  10.       double precision w(MAX_SIZE)
  11.       integer n, plusminus, inc
  12.       
  13.       double precision u, v, tpi
  14.       TPI = 2 * 3.14159265358979323846
  15.  
  16.       if( inc .ne . 1) goto 11
  17.       if ( plusminus .eq. -1) then
  18.     l = (n+2)/2
  19. c$doacross local(k,u,v), share(w)
  20.     do k = 1, l
  21.         u = 0.0
  22.         v = 0.0
  23.         do i = 1, n
  24.         u = u + a(i) * cos((k-1) * (i-1) * tpi / float(n))
  25.         v = v - a(i) * sin((k-1) * (i-1) * tpi / float(n))
  26.         end do
  27.         w((2*k-2)+1) = u
  28.         w((2*k-1)+1) = v
  29.     end do
  30. c$doacross
  31.     do i = 1, 2*l
  32.         a(i) = w(i)
  33.     end do
  34.       elseif (plusminus .eq. 1) then
  35.     print *,'NOT Implemented'
  36.       else
  37.     print *, 'Error first argument of DFT1DU should -1 or +1'
  38.     stop
  39.  11    print *, 'This Simple version support only INC = 1'
  40.     stop
  41.       endif
  42.  
  43.       return
  44.       end
  45.  
  46.  
  47.       subroutine zft1d( plusminus, n, a, inc)
  48.       implicit double precision ( a-h, o-z)
  49.       double complex a(*)
  50.       double complex w(MAX_SIZE)
  51.       integer n, plusminus
  52.       
  53.       double complex cp
  54.       double precision re, im, tpi
  55.  
  56.       tpi = 2.*3.14159265358979323846
  57.  
  58.     do i = 1, n
  59.         w(i) = 0.0
  60.         do k = 1, n
  61.         im = mod(((i-1)*(k-1)), n)
  62.         im = dfloat( plusminus*im)*tpi / dfloat(n) 
  63.         re = cos( im)
  64.         im = sin( im)
  65.         cp = dcmplx( re, im)
  66.         w(i) = w(i) + cp*a(inc*(k-1)+1)
  67.         end do
  68.     end do
  69.     do i = 1, n
  70.         a(inc*(i-1)+1) = w(i)
  71.     end do
  72.  
  73.       return
  74.       end
  75.  
  76.  
  77.  
  78. *****************************************************************************
  79. *                                        *
  80. *    2D (slow) Fourier Transform                        *
  81. *                                        *
  82. *****************************************************************************
  83.       subroutine zft2d( plusminus, n1, n2, w, ldw)
  84.       double complex w(ldw,n2)
  85.       double complex save( MAX_SIZE)
  86.       integer n1, n2, lda, ldw, plusminus
  87.  
  88.       call zfft1di( n1, save)
  89. c$doacross local(j)
  90.       do j = 1, n2
  91.     call zfft1d( plusminus, n1, w(1,j), 1, save)
  92.       end do
  93.  
  94.       call zfft1di( n2, save)
  95. c$doacross local(i,j,k,re,im,cp)
  96.       do i = 1, n1
  97.     call zfft1d( plusminus, n2, w(i,1), ldw, save)
  98.       end do
  99.  
  100.       return
  101.       end
  102.  
  103.       subroutine dft2du( plusminus, n1, n2, w, ldw)
  104.       Double Precision w(ldw,n2)
  105.       Double Precision save( 2*MAX_SIZE), cp(2*MAX_SIZE,MAX_THREADS)
  106.       integer n1, n2, lda, ldw, plusminus, l
  107.  
  108.       l = ((n1+2)/2)*2
  109.       if (ldw .lt. l) then
  110.     PRINT *, 'ERROR Leading Dimension is not sufficient',short(l),short(ldw)
  111.       endif
  112.       if (plusminus .eq. -1) then
  113.     call dfft1di( n1, save)
  114. c$doacross local(j)
  115.     do j = 1, n2
  116.         call dfft1du( plusminus, n1, w(1,j), 1, save)
  117.     end do
  118.  
  119.     call zfft1di( n2, save)
  120. c$doacross local(i,j,my), share(cp)
  121.     do i = 1, l, 2
  122.         my = 1
  123. c$        my = mp_my_threadnum()+1
  124.         do j = 1, n2
  125.         cp(2*j-1,my) = w(i+0,j)
  126.         cp(2*j-0,my) = w(i+1,j)
  127.         end do
  128.         call zfft1d( plusminus, n2, cp(1,my), 1, save)
  129.         do j = 1, n2
  130.         w(i+0,j) = cp(2*j-1,my) 
  131.         w(i+1,j) = cp(2*j-0,my)
  132.         end do
  133.     end do
  134.       else
  135.     call zfft1di( n2, save)
  136. c$doacross local(i,j,my)
  137.     do i = 1, l, 2
  138.         my = 1
  139. c$        my = mp_my_threadnum()+1
  140.         do j = 1, n2
  141.         cp(2*j-1,my) = w(i+0,j)
  142.         cp(2*j-0,my) = w(i+1,j)
  143.         end do
  144.         call zfft1d( plusminus, n2, cp(1,my), 1, save)
  145.         do j = 1, n2
  146.         w(i+0,j) = cp(2*j-1,my) 
  147.         w(i+1,j) = cp(2*j-0,my)
  148.         end do
  149.     end do
  150.  
  151.     call dfft1di( n1, save)
  152. c$doacross local(j)
  153.     do j = 1, n2
  154.         call dfft1du( plusminus, n1, w(1,j), 1, save)
  155.     end do
  156.       end if
  157.       return
  158.       end
  159.       
  160.  
  161. *****************************************************************************
  162. *
  163. *    3D (slow) Fourier Transform
  164. *
  165. *****************************************************************************
  166.  
  167.       subroutine zft3d(plusminus,n1,n2,n3,w,ld1,ld2)
  168.       double complex w(ld1,ld2,n3)
  169.       double complex save( MAX_SIZE)
  170.       integer n1,n2,n3,ld1,ld2,plusminus
  171. c
  172. c   transform along X  first ...
  173. c
  174.       call zfft1di( n1, save)
  175. c$doacross local(i,j,k)
  176.       do k = 1, n3
  177.        do j = 1, n2
  178.         call zfft1d( plusminus, n1, w(1,j,k), 1, save)
  179.        end do
  180.       end do
  181. c
  182. c   transform along Y then ...
  183. c
  184.       call zfft1di( n2, save)
  185. c$doacross local(i,j,k)
  186.       do k = 1,n3
  187.        do i = 1,n1
  188.         call zfft1d( plusminus, n2, w(i,1,k), ld1, save)
  189.        end do
  190.       end do
  191. c
  192. c   transform along Z finally ...
  193. c
  194.       call zfft1di( n3, save)
  195. c$doacross local(i,j,k)
  196.       do i = 1, n1
  197.        do j = 1, n2
  198.         call zfft1d( plusminus, n3, w(i,j,1), ld1*ld2, save)
  199.        end do
  200.       end do
  201.  
  202.       return
  203.       end
  204.  
  205.       subroutine dft3du(plusminus,n1,n2,n3,w,ld1,ld2)
  206.       Double Precision w(ld1,ld2,n3)
  207.       Double Precision save(2*MAX_SIZE), cp(MAX_SIZE,MAX_THREADS)
  208.       integer n1,n2,n3,ld1,ld2,plusminus, l
  209.  
  210.       if( plusminus .ne. -1) then
  211.     print *, 'Option non implemented YET'
  212.       end if
  213.       l = ((n1+2)/2)*2
  214.       if( ld1 .lt. l ) then
  215.     print *, 'Leading dimension not sufficient'
  216.       endif
  217. c
  218. c   trandform along X  first ...
  219. c
  220.       call dfft1di( n1, save)
  221. c$doacross local(i,j,k)
  222.       do k = 1, n3
  223.        do j = 1, n2
  224.         call dfft1du( plusminus, n1, w(1,j,k), 1, save)
  225.        end do
  226.       end do
  227. c
  228. c   trandform along Y then ...
  229. c
  230.       call zfft1di( n2, save)
  231. c$doacross local(i,j,k,my), share(cp)
  232.       do k = 1,n3
  233.     my = 1
  234. c$    my = mp_my_threadnum()+1
  235.     do i = 1, l, 2
  236.         do j = 1, n2
  237.         cp(2*j-1,my) = w(i+0,j,k)
  238.         cp(2*j-0,my) = w(i+1,j,k)
  239.         end do
  240.         call zfft1d( plusminus, n2, cp(1,my), 1, save)
  241.         do j = 1, n2
  242.         w(i+0,j,k) = cp(2*j-1,my) 
  243.         w(i+1,j,k) = cp(2*j-0,my)
  244.         end do
  245.     end do
  246.       end do
  247. c
  248. c   trandform along Z finally ...
  249. c
  250.       call zfft1di( n3, save)
  251. c$doacross local(i,j,k,my)
  252.       do j = 1,n2
  253.     my = 1
  254. c$    my = mp_my_threadnum()+1
  255.     do i = 1, l, 2
  256.         do k = 1, n3
  257.         cp(2*k-1,my) = w(i+0,j,k)
  258.         cp(2*k-0,my) = w(i+1,j,k)
  259.         end do
  260.         call zfft1d( plusminus, n3, cp(1,my), 1, save)
  261.         do k = 1, n3
  262.         w(i+0,j,k) = cp(2*k-1,my)
  263.         w(i+1,j,k) = cp(2*k-0,my)
  264.         end do
  265.     end do
  266.       end do
  267.  
  268.       return
  269.       end
  270.  
  271. *****************************************************************************
  272. *                                        *
  273. *    1D (slow) Fourier Transform                        *
  274. *                                        *
  275. *****************************************************************************
  276.       subroutine sft1du( plusminus, n, a, inc)
  277.       real a(*)
  278.       real w(MAX_SIZE)
  279.       integer n, plusminus, inc
  280.       
  281.       real u, v, tpi
  282.       TPI = 2 * 3.14159265358979323846
  283.  
  284.       if( inc .ne . 1) goto 11
  285.       if ( plusminus .eq. -1) then
  286.     l = (n+2)/2
  287. c$doacross local(k,u,v), share(w)
  288.     do k = 1, l
  289.         u = 0.0
  290.         v = 0.0
  291.         do i = 1, n
  292.         u = u + a(i) * cos((k-1) * (i-1) * tpi / float(n))
  293.         v = v - a(i) * sin((k-1) * (i-1) * tpi / float(n))
  294.         end do
  295.         w((2*k-2)+1) = u
  296.         w((2*k-1)+1) = v
  297.     end do
  298. c$doacross
  299.     do i = 1, 2*l
  300.         a(i) = w(i)
  301.     end do
  302.       elseif (plusminus .eq. 1) then
  303.  11    print *,'NOT Implemented'
  304.       endif
  305.  
  306.       return
  307.       end
  308.  
  309.  
  310.       subroutine cft1d( plusminus, n, a, inc)
  311.       implicit real ( a-h, o-z)
  312.       complex a(*)
  313.       complex w(MAX_SIZE)
  314.       integer n, plusminus
  315.       
  316.       complex cp
  317.       real re, im, tpi
  318.  
  319.       tpi = 2.*3.14159265358979323846
  320.  
  321.     do i = 1, n
  322.         w(i) = 0.0
  323.         do k = 1, n
  324.         im = mod(((i-1)*(k-1)), n)
  325.         im = float( plusminus*im)*tpi / float(n) 
  326.         re = cos( im)
  327.         im = sin( im)
  328.         cp = cmplx( re, im)
  329.         w(i) = w(i) + cp*a(inc*(k-1)+1)
  330.         end do
  331.     end do
  332.     do i = 1, n
  333.         a(inc*(i-1)+1) = w(i)
  334.     end do
  335.  
  336.       return
  337.       end
  338.  
  339.  
  340.  
  341. *****************************************************************************
  342. *                                        *
  343. *    2D (slow) Fourier Transform                        *
  344. *                                        *
  345. *****************************************************************************
  346.       subroutine cft2d( plusminus, n1, n2, w, ldw)
  347.       complex w(ldw,n2)
  348.       complex save( MAX_SIZE)
  349.       integer n1, n2, lda, ldw, plusminus
  350.  
  351.       call cfft1di( n1, save)
  352. c$doacross local(j)
  353.       do j = 1, n2
  354.     call cfft1d( plusminus, n1, w(1,j), 1, save)
  355.       end do
  356.  
  357.       call cfft1di( n2, save)
  358. c$doacross local(i,j,k,re,im,cp)
  359.       do i = 1, n1
  360.     call cfft1d( plusminus, n2, w(i,1), ldw, save)
  361.       end do
  362.  
  363.       return
  364.       end
  365.  
  366.       subroutine sft2du( plusminus, n1, n2, w, ldw)
  367.       real w(ldw,n2)
  368.       real save( 2*MAX_SIZE), cp(2*MAX_SIZE,MAX_THREADS)
  369.       integer n1, n2, lda, ldw, plusminus, l
  370.  
  371.       l = ((n1+2)/2)*2
  372.       if (ldw .lt. l) then
  373.     PRINT *, 'ERROR Leading Dimension is not sufficient',short(l),short(ldw)
  374.       endif
  375.       if (plusminus .eq. -1) then
  376.     call sfft1di( n1, save)
  377. c$doacross local(j)
  378.     do j = 1, n2
  379.         call sfft1du( plusminus, n1, w(1,j), 1, save)
  380.     end do
  381.  
  382.     call cfft1di( n2, save)
  383. c$doacross local(i,j,my), share(cp)
  384.     do i = 1, l, 2
  385.         my = 1
  386. c$        my = mp_my_threadnum()+1
  387.         do j = 1, n2
  388.         cp(2*j-1,my) = w(i+0,j)
  389.         cp(2*j-0,my) = w(i+1,j)
  390.         end do
  391.         call cfft1d( plusminus, n2, cp(1,my), 1, save)
  392.         do j = 1, n2
  393.         w(i+0,j) = cp(2*j-1,my) 
  394.         w(i+1,j) = cp(2*j-0,my)
  395.         end do
  396.     end do
  397.       else
  398.     call cfft1di( n2, save)
  399. c$doacross local(i,j,my)
  400.     do i = 1, l, 2
  401.         my = 1
  402. c$        my = mp_my_threadnum()+1
  403.         do j = 1, n2
  404.         cp(2*j-1,my) = w(i+0,j)
  405.         cp(2*j-0,my) = w(i+1,j)
  406.         end do
  407.         call cfft1d( plusminus, n2, cp(1,my), 1, save)
  408.         do j = 1, n2
  409.         w(i+0,j) = cp(2*j-1,my) 
  410.         w(i+1,j) = cp(2*j-0,my)
  411.         end do
  412.     end do
  413.  
  414.     call sfft1di( n1, save)
  415. c$doacross local(j)
  416.     do j = 1, n2
  417.         call sfft1du( plusminus, n1, w(1,j), 1, save)
  418.     end do
  419.       end if
  420.       return
  421.       end
  422.       
  423.  
  424. *****************************************************************************
  425. *
  426. *    3D (slow) Fourier Transform
  427. *
  428. *****************************************************************************
  429.  
  430.       subroutine cft3d(plusminus,n1,n2,n3,w,ld1,ld2)
  431.       complex w(ld1,ld2,n3)
  432.       complex save( MAX_SIZE)
  433.       integer n1,n2,n3,ld1,ld2,plusminus
  434. c
  435. c   transform along X  first ...
  436. c
  437.       call cfft1di( n1, save)
  438. c$doacross local(i,j,k)
  439.       do k = 1, n3
  440.        do j = 1, n2
  441.         call cfft1d( plusminus, n1, w(1,j,k), 1, save)
  442.        end do
  443.       end do
  444. c
  445. c   transform along Y then ...
  446. c
  447.       call cfft1di( n2, save)
  448. c$doacross local(i,j,k)
  449.       do k = 1,n3
  450.        do i = 1,n1
  451.         call cfft1d( plusminus, n2, w(i,1,k), ld1, save)
  452.        end do
  453.       end do
  454. c
  455. c   transform along Z finally ...
  456. c
  457.       call cfft1di( n3, save)
  458. c$doacross local(i,j,k)
  459.       do i = 1, n1
  460.        do j = 1, n2
  461.         call cfft1d( plusminus, n3, w(i,j,1), ld1*ld2, save)
  462.        end do
  463.       end do
  464.  
  465.       return
  466.       end
  467.  
  468.       subroutine sft3du(plusminus,n1,n2,n3,w,ld1,ld2)
  469.       real w(ld1,ld2,n3)
  470.       real save(2*MAX_SIZE), cp(MAX_SIZE,MAX_THREADS)
  471.       integer n1,n2,n3,ld1,ld2,plusminus, l
  472.  
  473.       if( plusminus .ne. -1) then
  474.     print *, 'Option non implemented YET'
  475.       end if
  476.       l = ((n1+2)/2)*2
  477.       if( ld1 .lt. l ) then
  478.     print *, 'Leading dimension not sufficient'
  479.       endif
  480. c
  481. c   transform along X  first ...
  482. c
  483.       call sfft1di( n1, save)
  484. c$doacross local(i,j,k)
  485.       do k = 1, n3
  486.        do j = 1, n2
  487.         call sfft1du( plusminus, n1, w(1,j,k), 1, save)
  488.        end do
  489.       end do
  490. c
  491. c   transform along Y then ...
  492. c
  493.       call cfft1di( n2, save)
  494. c$doacross local(i,j,k,my), share(cp)
  495.       do k = 1,n3
  496.     my = 1
  497. c$    my = mp_my_threadnum()+1
  498.     do i = 1, l, 2
  499.         do j = 1, n2
  500.         cp(2*j-1,my) = w(i+0,j,k)
  501.         cp(2*j-0,my) = w(i+1,j,k)
  502.         end do
  503.         call cfft1d( plusminus, n2, cp(1,my), 1, save)
  504.         do j = 1, n2
  505.         w(i+0,j,k) = cp(2*j-1,my) 
  506.         w(i+1,j,k) = cp(2*j-0,my)
  507.         end do
  508.     end do
  509.       end do
  510. c
  511. c   transform along Z finally ...
  512. c
  513.       call cfft1di( n3, save)
  514. c$doacross local(i,j,k,my)
  515.       do j = 1,n2
  516.     my = 1
  517. c$    my = mp_my_threadnum()+1
  518.     do i = 1, l, 2
  519.         do k = 1, n3
  520.         cp(2*k-1,my) = w(i+0,j,k)
  521.         cp(2*k-0,my) = w(i+1,j,k)
  522.         end do
  523.         call cfft1d( plusminus, n3, cp(1,my), 1, save)
  524.         do k = 1, n3
  525.         w(i+0,j,k) = cp(2*k-1,my)
  526.         w(i+1,j,k) = cp(2*k-0,my)
  527.         end do
  528.     end do
  529.       end do
  530.  
  531.       return
  532.       end
  533.  
  534.